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We introduce a multi-discontinuity algorithm for the efficient global update of world-line configura- 
tions in Monte Carlo simulations of interacting quantum systems. This algorithm is a generalization 
of the two-discontinuity algorithms introduced in Refs. [N. Prokof'ev, B. Svistunov, and I. Tupit- 
syn, Phys. Lett. A 238, 253 (1998)] and [O. Syljuasen and A. Sandvik, Phys. Rev. E 66, 046701 
(2002)] . This generalization is particularly effective for studying Bose-Einstein condensates (BEC) 
of composite particles. In particular, we demonstrate the utility of the generalized algorithm by 
simulating a Hamiltonian for an S — 1 anti-ferromagnet with strong uniaxial single-ion anisotropy. 
The multi-discontinuity algorithm not only solves the freezing problem that arises in this limit, 
but also allows the efficient computing of the off-diagonal correlator that characterizes a BEC of 
composite particles. 

PACS numbers: 02.70.Ss,75.30.Kz 



Introduction: Cooperative phenomena in quantum sys- 
tems (e.g., superfluidity, Bose-Einstein condensation, su- 
perconductivity, magnetism, etc.) have fascinated physi- 
cists for decades. The study of these phenomena often 
requires one to consider interaction regimes for which tra- 
ditional perturbative approaches do not work. World-line 
quantum Monte-Carlo (QMC) simulations based on the 
Feynman path integral formulation have demonstrated 
to be a very powerful method for studying systems in 
these regimes [TJ [2] . The main reasons are that the QMC 
method is unbiased, and its results are exact within sta- 
tistical error. Therefore, it is imperative to expand the 
range of applicability of this method by developing QMC 
algorithms for problems that can not be simulated with 
state-of-the-art numerical techniques. 

Although the worm algorithm [3] and the directed-loop 
algorithm (DLA) [3] , are the most efficient strategies for 
global updates in the QMC simulations, they have diffi- 
culties when applied to systems exhibiting Bose-Einstein 
condensates (BECs) of composite particles. BECs of 
composite particles appear in different contexts, such as 
bosonic gases of two species (a and b) of ultra cold atoms 
[SJ (H] and homonuclear bosonic gases ( 23 Na2 [7], 133 Cs2 
[5], etc.). Another simple example of BECs of composite 
particles is the ferro-nematic magnetic ordering that ap- 
pears in quantum magnets with single-ion anisotropy [9] . 
In general, high spin systems (S > 1) whose total mag- 
netization along a particular axis is conserved (it will be 
our z-axis) can exhibit a local order parameter of the 
form ((S + ) n ) (2 < n < 25) in the absence of magnetic 
and any lower moments, i.e., ((S + ) m ) = for m < n. 
These order parameters correspond to different moments 
of a distribution of electric currents. For instance, n = 2 
corresponds to quadrupolar or nematic order, and n = 3 
corresponds to octupolar order. After mapping the spin 
system into a gas of bosons, the phase with order pa- 
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rameter ((S + ) n ) ^ is mapped into a BEC of composite 
particles containing rt-bosons. For instance, nematic or- 
dering is mapped into a BEC of pairs of bosons [5] ■ 

QMC simulations sample and update the d + 1 dimen- 
sional world-line configurations of the quantum system 
under consideration (d is the spatial dimension). Stan- 
dard algorithms insert a pair of discontinuities head and 
tail in the world-line configuration. The head propagates 
until it meets the tail [2]. The head changes the local 
states of the world-line configuration along its way. As a 
result, states on the loop drawn by the head's trajectory 
are globally updated. For bosonic systems, the disconti- 
nuities typically correspond to the creation, at , or annihi- 
lation, a, of a particle. For QMC simulations of mixtures 
of bosonic gases [5J HDHTS] . including a-type and b- 
type bosons, the discontinuities can be of the form atb 
or b^a that exchange the type of boson [12]. However, 
we will see that two discontinuities are not enough to 
perform efficient simulations in certain cases because of 
the so-called "freezing problem" . To solve this problem, 
algorithms with two or more pairs of discontinuities cor- 
responding to a, a\b or have been considered as a gen- 
eralization of the worm algorithm (e.g. Refs. [51 [T3] ) . 
Furthermore the two-discontinuity algorithms often suf- 
fer from a serious slowing down problem for calculations 
of off-diagonal correlators like (aJ.(r)6 T .(r)&J,,(r / )a r .'(T / )), 
which characterize the condensation of composite parti- 
cles. 

In this paper, we present a method called the multi- 
discontinuity algorithm (MDA) that incorporates an ar- 
bitrary number of discontinuities in order to solve the 
freezing problem. We will show that the MDA can be 
used for simulating BECs of composite particles that can- 
not be simulated efficiently with two-discontinuity algo- 
rithms. Although the MDA can be applied to many dif- 
ferent situations, we will consider a simple example of 
an S = 1 antiferromagnetic (J > 0) Hciscnberg model 
with external magnetic field, B, and single-ion uniaxial 
anisotropy D. We will see that the MDA successfully 
explores the ferronematic phase of this model, while the 
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DLA fails because of the freezing problem. Since one ad- 
ditional discontinuity is enough for simulating the S = 1 
system, we restrict to three discontinuities in this paper. 
The extension to the algorithm to more than three dis- 
continuities is straightforward. 

Model: We consider the spin-one Hamiltonian 



u 



J S r - S r , +Dj2Sf -Bj2 S r, (1) 



(r ,r') 



where the summation of (r, r') runs over all the 
nearest-neighbor bonds and the model is defined on 
d-dimensional hypercubic lattices [TBJ. We restrict 
the model to the case of easy-axis single-ion anisotropy 
D < 0. In the absence of a magnetic field (B = 0), the 
ground state of % exhibits Ising-like antiferromagnetic 
ordering along the z axis (AFM-z). This AFM-z phase 
is destroyed by the application of a strong enough mag- 
netic field B via a first-order quantum phase transition to 
a new phase that exhibits ferronematic (FN) ordering [9] . 
A further increase of the magnetic field leads to a second 
order phase transition to the fully polarized phase. As 
we will see later, the MDA is necessary for finding and 
studying the FN phase that is characterized by the FN 
susceptibility, 



Xfn 
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Algorithm: World-line configurations in the d + 1 di- 
mensional space introduced in the Feynman path integral 
formulation are sampled by the QMC method [2J. Here 
we explain the MDA as a generalization of the DLA by 
using model ([I]) as an example of its application. When 
the DLA is simply applied to model Q for J > 0, it 
suffers from the negative sign problem. To avoid the 
negative sign problem, we introduce a spin rotation by 
7r along the z axis (S^.' y —> —S*' v ) for the B sublattice. 
The expression of the Hamiltonian in the new basis is 

T~L = J ^ (S r S r r -f~ S y S^,/ S r S r r) 

(r,r') 

+DY J S?-BY J S z r . (2) 

r r 

We will start by applying the DLA to see how the freez- 
ing problem arises for \D\ ^> J and B > J. Figure 
[lja) shows a typical world-line configuration appearing 
at \D\ ^> J and low temperature. Because of the strong 
anisotropy, \D\ 3> J, the states \S Z — ±1) are the vast 
majority, while the \S Z =0) states rarely appear in the 
world- line configurations. The density of these minor- 
ity states is ~ J 2 /D (note that the volume of the d + 1 
dimensional space has units of inverse energy so the den- 
sity of states has units of energy). Moreover, the \S Z = 0) 
states typically appear forming bubbles that correspond 
to quantum fluctuations. These quantum fluctuations 
are crucial for the stabilization of the FN ordering [3]. 



We will see that the bubbles of IS" 2 = 0) states are the 
main source of difficulty for the DLA. 

In the DLA algorithm, the world-line configurations 
are globally updated by introducing a pair of disconti- 
nuities called worm and their scattering objects called 
vertices. The world-line configurations are updated ac- 
cording to the following protocol: 

Step 1 Place vertices stochastically (purple double line 
in Fig. [ijb)). The density of vertices at ((r, r'),r) 
is a function of (S z (t), S z ,(t)). 

Step 2 Choose a point using a uniform random distri- 
bution, create a pair of discontinuities (worm), and 
assign one discontinuity to be the head (triangle 
object in Fig. |TJb) ) and the other to be the tail 
(circle object in Fig. [ljb)). 

Step 3 Move the head until it hits a vertex or the tail. 
If it hits the tail, go to Step 3-1, otherwise, Step 
3-2. 

Step 3-1 Let the worm annihilate, and go back to Step 
2, or go to Step 4. 

Step 3-2 Choose its new direction stochastically, and go 
to Step 3. 

Step 4 Erase all vertices after creating and annihilating 
the worm several times. 

More details of the DLA can be found in Ref. [H0]. The 
worm is generated by a "source term", — rjoQ, that is 
added to the Hamiltonian (770 > 0) . For the case of our 
model ([2]), we choose Q as 



Q = 



{St 



f + (^-) 2 
2 



(3) 



The discontinuity is created by any of the four terms S£ , 
S^, (S^) 2 , or (5~) 2 . We sample world-line configura- 
tions whose Boltzmann weights are 0(t]q) and 0(r^). In 
other words, the allowed world-line configurations have 
only one warm or no discontinuity at all. Because of 
the external field B > J, the head created by (<S* + ) 2 is 
trapped locally in a bubble, as shown in Fig. [2j The 
probability of escaping from such a trap is exponen- 
tially small in the number of vertices corresponding the 
B term on the way to escape, and it becomes negligi- 
bly small because the typical number of the vertices is 
0(DB/ J 2 ) » O(l). As a result, it takes an exponen- 
tially long time to complete a cycle of creation and an- 
nihilation of a worm. This trapping problem makes the 
DLA inefficient for simulating the FN phase of H. 

The basic difficulty of the the DLA can be eliminated 
by considering an enlarged Hilbert space via the inclu- 
sion of a third discontinuity. We add a new source term, 
—rjiQ to % — rjoQ, to generate a third discontinuity that 
has a weight 771, which can be different from 770 , and 
sample the world-line configurations whose weights are 
0(VoVi)i 0(i] 2 rii), or 0(ry 2 ?yj). Accordingly, we consider 
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FIG. 1: (Color online) (a) Typical world-line configuration 
for a chain of six sites (d = 1), a large \D\ 3> J ratio, and 
very low temperature j3 J 3> 1. Red, gray, and blue solid lines 
correspond to S z — —1, 0, and +1, respectively. The stars (*) 
mark all the bubbles that represent quantum fluctuations, 
(b) A world-line configuration with a pair of discontinuities 
called a worm in the directed-loop algorithm. The purple 
horizontal double lines are the scattering objects (vertices) 
for the discontinuity. 
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FIG. 2: (Color online) Local trapping of an (S + ) 2 discon- 
tinuity at a bubble. The vertices are not illustrated here for 
simplicity. 



three new processes that will be called "fusion" , "fission" 
and "swap" (see Fig. [|). These new processes are ex- 
ecuted when the head stops at a newly introduced one- 
body vertex. Although this new vertex is supposed to be 
created with a uniform distribution, in practice we cre- 
ate a virtual new vertex [T71 [T5] to avoid the unnecessary 
proliferation of vertices and the introduction of new free 
parameters. Fission corresponds to the process by which 
one discontinuity decays into two discontinuities. We rep- 
resent these processes by indicating the source operators 



that create the initial and the final states: [(S + ) 2 — > S + 
and S+], [(S*-) 2 S- and S~], [S+ -> (S+) 2 and S~], 
and [S~ — > (S* - ) 2 and S + ). The inverse processes cor- 
respond to fusions because they combine two disconti- 
nuities into a single one. The swap corresponds to the 
process of exchanging the head and one of the two tails. 
To include these three new processes, we need to replace 
Steps 3, 3-1, and 3-2 of the conventional DLA by the 
following steps: 

Step 3 Calculate I = — ln(i?)/(2^i). R is a random 
number with uniform distribution in the interval 
(0,1]. (See the Appendix for details of virtual 
placement of vertices.) If Z > tq, where tq is the 
time interval between the head and the next ob- 
ject the head finds on its way (a vertex or another 
discontinuity), go to Step 3-1, otherwise go to Step 
3-2. 

Step 3-1 Let the head go to the next object. If the next 
object is a vertex, go to Step 3-1-1, otherwise go to 
Step 3-1-2. 

Step 3-1-1 Choose a new direction of the head stochas- 
tically, and go back to Step 3. 

Step 3-1-2 If the total number of discontinuities is nr> = 
2, go to Step 3-1-2-1, otherwise go to Step 3-1-2-2. 

Step 3-1-2-1 Let the discontinuities annihilate each 
other, and go back to Step 2, or go to Step 4. 

Step 3-1-2-2 Let the two discontinuities fuse and go to 
Step 3. 

Step 3-2 Move the discontinuity through the time in- 
terval I, If ri£, = 2, go to Step 3-2-1, otherwise 
(n D =3) go to Step 3-2-2. 

Step 3-2-1 Let the head fission, and go to Step 3. 

Step 3-2-2 Let the head freeze, exchange discontinuities 
as a new head to move randomly (swap) , and go to 
Step 3. 

The value of the free parameter 771 is fixed so that the 
average path lengths of the head with and without the 
third discontinuity are roughly the same. 

Example of application: Figure [4] shows the FN sus- 
ceptibility xfn °f model ^ defined on a ring (d = 1) of 
L = 5 spins and for Hamiltonian parameters D/J = 0, 
and B/J = 0.3. The three curves are the results ob- 
tained by using MDA, DLA and exact diagonalization 
(ED). The agreement between MDA and ED results is 
excellent. However, there is a rather large deviation be- 
tween the DLA and ED results, although the agreement 
is very good for other physical quantities (e.g., energy 
and uniform magnetization) that are not shown. This 
deviation confirms that DLA does not obtain an accurate 
estimate of the off-diagonal pair-pair correlation function 

(s+ + (T + )S+ + (T + )S-_(T-)S-_(T-)y that determines the 
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FIG. 3: (Color online) Example of a sequence of processes 
involving "fusion" , "fission" , and "swap" , that allows the head 
to escape from the trap in the multi-discontinuity algorithm. 



FN susceptibility because of a serious slowing down prob- 
lem. This correlation function is obtained by counting 
the number of events in which only two discontinuities 
corresponding to (S + ) 2 and (S~) 2 exist in whole space, 
with the (S + ) 2 discontinuity located at (r,r) = (r + ,r + ) 
and the (S~) 2 discontinuity located at (r,r) = (r_,r_) 
[12] . As it is clear from Eq. xfn is the integration 
of the sum of this correlation function over all the possi- 
ble values of the coordinates (r + ,r + ) and (r_,r_). We 
sample the whole length of the trajectory of the (S l+ ) 2 
and (S~) 2 discontinuities traveling when no other discon- 
tinuity exists. As we explained in the previous section, 
the DLA suffers from a serious slowing down problem 
because the head of a worm made of a pair of (S + ) 2 
and (S~) 2 discontinuities cannot go through a bubble of 
S z = states. This limitation precludes the correct esti- 
mation of xfn, which is the crucial physical quantity for 
identifying the FN phase in a finite size system. 
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FIG. 4: (Color online) Ferronematic susceptibility, Xfn, of 
model |2| defined on a ring (d = 1) of L — 5 spins. The Hamil- 
tonian parameters are D/J = 0, and B/J = 0.3. The three 
curves are the results obtained with the multi-discontinuity al- 
gorithm (MDA), the original directed- loop algorithm (DLA) 
and exact diagonalization (ED). We performed 10 7 Monte- 
Carlo sweeps for both MDA and DLA. 



Figure [5] shows the temperature dependence of the con- 
veniently scaled FN susceptibility, xfn£ _2+ '' , for the 
case of a simple cubic lattice and Hamiltonian param- 
eters D/J = —8 and B/J = 6.25. The common crossing 
point among the curves obtained with different system 
sizes indicates that there is a finite temperature transi- 
tion between the low temperature FN phase and a high 
temperature paramagnetic phase. The order parameter 
of the FN phase is ((S^) 2 ) , which in bosonic language can 
be interpreted as a Bose-Einstein condensate of pairs of 
bosons. Thus, if continuous, the phase transition should 
belong to the 3D XY universality class. This is indeed 
the case according to the scaling analysis shown in Fig. [5] 
there is a common crossing point when we choose the 
value of r\ = 0.038 that corresponds to the 3D XY uni- 
versality class. 
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FIG. 5: (Color online) Finite-size scaling of ferronematic sus- 
ceptibility for H defined on a simple cubic lattice [d = 3) with 
periodic boundary conditions. The Hamiltonian parameters 
are at D/J = —8, and B/J — 6.25. The estimated criti- 
cal temperature is T c = 0.057. We use the critical exponent 
r\ = 0.038 for the 3D XY universality class [Tg]. 

Conclusions: In summary, we have introduced a MDA 
that overcomes some limitations of the two-discontinuity 
algorithms. The MDA helps to solve the freezing prob- 
lem that often arises in the simulation of strongly inter- 
acting quantum systems. The example that we discussed 
is just one out of many. Moreover the two-discontinuity 
algorithms do not satisfy ergodicity for a large class of 
physically relevant models. In other words, there are 
physically possible world-line configurations that have 
not been sampled. For instance, the MDA is also neces- 
sary for simulating any quantum spin system with cubic 
single-ion anisotropy [2j . This is a quite common case for 
several Mott insulating materials such as EuTiC>3 [20 . 

The extension of the three-discontinuity algorithm to 
Y-discontinuity is straightforward, and it is necessary 
for describing certain physical situations such as S > 2 
antiferromagnets with cubic single-ion anisotropy. The 
only qualitative difference for Y > 3 is the possibility 
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of creating or annihilating a new pair of discontinuities 
in world-line configurations that already contain a finite 
number of discontinuities. For efficient simulations with 
N > 3, we choose the swap or the fission stochastically 
when the head stops at a virtual vertex with ri£> < N. 
(We choose the fission with probability 1 when no = 2 
and N = 3.) Although we introduced the MDA as a 
generalization of the DLA in this paper, it is straightfor- 
ward to modify the MDA in such a way that it becomes 
a generalization of the worm algorithm. 

The successful application of the MDA to the S = 
1 Heiscnberg antiferromagnet with uniaxial single-ion 
anisotropy and Zeeman coupling to an external mag- 
netic field demonstrates the relevance of this method. We 
showed that the pair-condensation or ferronematic sus- 
ceptibility, xfn, obtained from simulations based on the 
MDA accurately reproduces the ED results. The same is 
not true for the DLA. Our MDA simulations also repro- 
duce the expected finite temperature second order phase 
transition between the paramagnetic and FN phases in a 
simple cubic lattice. 
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Appendix: Fission and fusion 

In this appendix, we consider the details of virtual 
placement of vertices for fission and fusion. We intro- 
duce a single-site vertex for the fission and fusion pro- 
cesses. The weight of the vertex W is uniform everywhere 
and independent from the local state of the world-line. 
We set W = 2r]i so that the acceptance probability of 
fusion is equal to unity. We start by considering the 
time-reversal symmetry of the algorithm that is a suffi- 



cient condition for satisfying detailed balance 2J. Here, 
we will only consider the time-reversal symmetry con- 
dition associated with the fission and fusion processes. 
We decompose the weights just before a fission or fusion 
operation as w a — w aP where w a is the weight of 
the initial configuration a. By using the partial weights, 
w a », we determine the transition probability from a to 
P as t a p = w a p I w a , so that the time-reversal symme- 
try condition, w a t a j3 = wptp a , is satisfied. Typically, we 
take w a p = wpa for simplicity. Figure [6] shows an ex- 
ample of w a p of the fission and fusion processes. The 
table is symmetric when W — 2rji and the time-reversal 
symmetry condition is satisfied. "Virtual placement of a 
vertex" means finding the next point along the world line 
where the head is scattered. The vertices are generated 
based on a Poisson distribution with density W. Namely, 
the time interval to the next vertex is — ln(i?)/VF, where 
R is random number between and 1: R €E (0, 1]. 
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FIG. 6: (Color online) Table of w a p related to the fission 
and fusion processes. Green crosses represent the newly intro- 
duced single-site vertices for the fission and fusion processes. 
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